Annealing a Magnetic Cactus into Phyllotaxis 
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The appearance of mathematical regularities in the disposition of leaves on a stem, scales on a 
pine-cone and spines on a cactus has puzzled scholars for millennia; similar so-called phyllotactic 
patterns are seen in self-organized growth, polypeptides, convection, magnetic flux lattices and ion 
beams. Levitov showed that a cylindrical lattice of repulsive particles can reproduce phyllotaxis 
under the (unproved) assumption that minimum of energy would be achieved by 2-D Bravais lat- 
tices. Here we provide experimental and numerical evidence that the Phyllotactic lattice is actually 
a ground state. When mechanically annealed, our experimental "magnetic cactus" precisely repro- 
duces botanical phyllotaxis, along with domain boundaries (called transitions in Botany) between 
different phyllotactic patterns. We employ a structural genetic algorithm to explore the more gen- 
eral axially unconstrained case, which reveals multijugate (multiple spirals) as well as monojugate 
(single spiral) phyllotaxis. 

PACS numbers: 87.10.-e, 68.65.-k, 89.75.Fb 



I. INTRODUCTION 

Symmetrical morphologies and regular patterns in liv- 
ing organisms (FigfTJ) have been credited with originating 
the idea of beauty, the notion of art as an imitation of 
nature, and humanity's first mathematical inquiries 
UJ. The fascinating symmetrical patterns of organs in 
plants, called phyllotaxis P~H2], were known to the Ro- 
mans (Pliny) and ancient Greeks (Theofrastus), while 
early recognitions are found in sources as ancient as the 
Text of the Pyramids [1] . Leonardo da Vinci [5] , Andrea 
Cesalpino, and Charles Bonnet [5] studied phyllotaxis in 
the modern era. Kepler proposed that the Fibonacci se- 
quence (1, 2, 3, 5, 8. . . ), where each term is the sum of 
the two preceeding ones [7|, describes these phyllotactic 
patterns. 

A discipline that thrived on multidisciplinary interac- 
tions 8J, phyllotaxis found its standard mathematical 
description when August and Louis Bravais [S] intro- 
duced the point lattice on a cylinder to represent the 
dispositions of leaves in 1837 (see Fig. [5]), thirteen years 
before August's seminal work on crystallography [10] . 
Unfortunately botanists neglected the work of the Bra- 
vais brothers, and it wasn't until Church rediscovered it 
eighty years later that more progress was achieved in the 
field [IT]. 

The geometrical description of cylindrical phyllotaxis 
relies, in the simplest case, on the phyllotactic lattice in- 
troduced by the Bravais brothers jTJ-3, 9J. It consists of 
a so-called generative spiral of divergence angle Q. We 
can visually decompose the resulting lattice in crossing 
spirals that join nearest neighbors, as in Fig. [2] which 
botanists call parastichies. It is a fundamental observa- 
tion (made first by Kepler) that the numbers n, m of 
crossing parastichies needed to cover the lattice are con- 
secutive terms of the standard Fibonacci sequence, or less 



frequently the variants obtained by changing the second 
term, also called Lucas numbers: 1, 3, 4, 7, 11. . . and 1, 4, 
5, 9. . . often referred to as second and third phyllotaxis. 
From that, one can prove that the divergence angle of the 
generative spirals in plants assumes values close to [3] [T2"] 



where p — 1, 2, 3 denotes first, second or third phyllotaxis 
and r = (l + \/5) /2 is the golden ratio. For more than 
one generative spiral ("multijugate" phyllotaxis), paras- 
tichies share a common divisor (n,m) — (kn' ,km'), k 
being the number of generative spirals [5] [T2] . Not un- 
like domain boundaries in crystals, plants show kinks be- 
tween domains, called transitions by botanists P~l . 

In the last 50 years, phyllotactic patterns have 
been seen or predicted outside of botany: polypep- 
tide chains [21 [15], tubular packings of spheres [16] . 
convection cells [17], layered superconductors [TS], self- 
assembled microstructures [Jj5], and cooled particle 
beams [201 HI]- While it is still debated whether such 
systems might shed light on botanical phyllotaxis, the 
occurrence of such mathematical regularities outside of 
botany is fascinating and leads to generalizations that - 
unlike quasistatic botany - allows for dynamics [22] . 

In a groundbreaking work Levitov recognized phyl- 
lotaxis in vortices of layered superconductors [18) . He 
next described how phyllotactic patterns represent states 
of minimal energy of a cylindrical lattice (that is of a lat- 
tice with cylindrical boundary conditions) of mutually 
repelling objects, the repulsion mimicking the interac- 
tions between spines, leaves, or seeds in plant morphol- 
ogy [23l [24]. Yet such a constraint to a lattice is ab- 
sent both in botany and in the physical systems to which 
this energetic model might apply, such as adatoms or 
low-density electrons on nanotubes and ions or dipolar 
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FIG. 1: Natural and Magnetic Cacti. A specimen of Mam- 
millaria elongata displaying a helical morphology ubiquitous 
to nature, a magnetic cactus of dipoles on stacked bearings, 
and a schematic of a wrapped Bravais lattice showing the an- 
gular offset (divergence angle) 0, and the axial separation a 
between particles. 



molecules in cylindrical traps. 

Following up on earlier work that focused on the dy- 
namics of rotons and solitons in physical phyllotactic 
systems (22] . we provide here a detailed experimental 
and numerical demonstration that Levitov's constraint 
is not necessary, and that the lowest energy states of 
repulsive particles in cylindrical geometries are indeed 
phyllotactic lattices. In addition, we describe the exper- 
imental and numerical generation of multijugate phyl- 
lotaxis, static kink-like domain boundaries between dif- 
ferent phyllotactic lattices, and unusual disordered yet 
reflection-symmetric structures that may be a static relic 
of soliton propagation. 

We show that when a "magnetic cactus" of magnets 
(spines) equally spaced on co-axial bearings (stem) with 
south poles all pointing outward is annealed, it precisely 
reproduces botanical phyllotaxis. When studied numeri- 
cally via a structural genetic algorithm, the fully uncon- 
strained case reveals both multijugate and mono jugate 
phyllotaxis. In addition to our macro-scale implementa- 
tion, such systems could also be created at the quantum 
level in nanotubes or cold atomic gases. 

In section II we describe the statics of repulsive par- 
ticles in cylindrical geometries. In section III we detail 
the experiment on the magnetic cactus. In Section IV we 
discuss the more general case of multijugate phyllotaxis. 

II. PHYLLOTAXIS OF REPULSIVE PARTICLES 
IN CYLINDRICAL GEOMETRIES 

In this section we will recall Levitov's model [331 121] 
and some of our own findings [32]. Following Levitov, 
let us assume that the lowest energy configuration for a 
set of particles with repulsive interactions, confined to a 
cylindrical shell of radius R, is a helix with a fixed angu- 
lar offset Cl between consecutive particles and a uniform 
axial spacing a, as in Fig. [I] (this so far unproved ansatz 
will be investigated later both numerically and experi- 
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FIG. 2: The Bravais lattice with cylindrical boundary con- 
ditions that defines a phyllotactic spiral. The cylinder axis 
is vertical, while the horizontal direction contains three cir- 
cumferential repeats. The solid line is the generative spiral: 
this one-dimensional Bravais lattice generates the full struc- 
ture. The dashed lines are the so-called parastichies or visible 
secondary spirals: they connect nearest neighbors on the sur- 
face of the cylinder. Adapted from A. Bravais and L. Bravais, 
1837 [9]. 



mentally). For a generic pair- wise repulsive interaction 
Vij between particles i and j, the energy of the helix is 

V = \ v i,j- Since the lattice structure is defined by 
Cl, we can write V{Cl), 

In Fig. |3] we plot V {CI) for various values of the ratio 
a/R: for specificity we employed a dipole dipole interac- 
tion = Pi-Pj/rfj - 3{pi ■ r it j){pj ■ rij)/r^, repulsive 
at the densities considered here. However, the following 
considerations only depend upon geometry and therefore 
apply to a vast range of reasonably behaved, long range 
repulsive interactions. 

When a/R ^> 1, the angle Cl — ir maximizes distance 
between neighboring particles and therefore V (Cl) has a 
minimum in ir. The angle between second nearest neigh- 
bors along the helix is 2n, which means that they face 
each other. And thus, as the density increases, inter- 
action between the facing second nearest neighbors be- 
comes predominant, and Cl = tt is not a minimum for 

V (Cl) anymore. If whe shift the helical angle from w, the 
repulsive interaction between second nearest neighbors is 
reduced, with minimal penalty from nearest neighbors. 
In terms of V {Cl), that means a local maximum Cl = n. 

This argument can be iterated for every commensu- 
rate winding that allows particles separated by j neigh- 
bors to face each other. As density increases further, 
the angles 27r/3 and 4tt/3 also become unfavorable due 
to third-neighbor interactions. Any commensurate spiral 
of divergence angle Cl = 2iri/j with i,j relatively prime 
corresponds to a configuration where each particle faces 
each jth neighbor. For every j there will be a value of 
a/R low enough such that Cl — 2iti/ j is a local maximum, 
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which we call a peak of rank j. 

The proliferation of peaks for increasing linear density 
is shown in Fig [3] We can see that at high density, peaks 
of equal rank are nearly degenerate; that is natural, since 
their principal defining energetic contribution arises from 
particles facing each other at a distance ja. The minima 
also become more nearly degenerate as the density in- 
creases. That can be explained intuitively, since for an- 
gles incommensurate to tt each particle "sees" the others 
as incommensurately smeared around the cylinder, and 
is therefore embedded in a nearly uniform background 
charge from the other particles. The degenerate energy 
of the ground state can be well approximated by an uni- 
form continuum distribution eo, whereas the energy of a 
peak of order j will be V (2m/ j) ~ v(ja) + eo, where v(r) 
is the energy of two particles facing at a distance r: for 
our dipole interaction v(ja) = p 2 /a 3 j 3 . 

The first step to calculate the degeneracy of our sys- 
tem at a given density, is to find the corresponding max- 
imum rank of the peaks. As all of the peaks of the same 
rank have the same energy, and appear in the spectrum 
together, we can focus on the emergence of 2ir/J. For 
a/R <C 1, this new peak will emerge when the distance 
between particles separated by a distance Ja equals that 
of particles separated by 2ttR/J. Therefore one finds for 
the maximum rank 

rr /27ri?n 

which as expected only depends on purely geometrical 
parameters. A little more tricky is to compute the de- 
generacy, given J. The set of all the peaks has the car- 
dinality of the class of all the fractions i/j, with i, j 
coprime and j < J. This can be considered as the dis- 
jointed union of other classes, called Farey classes of order 
j, defined as follows: Pj = {Cl = 2ni/j | for i, j coprime 
and i < j}, i.e. all fractions in lowest terms between 
and 1 whose denominators do not exceed j [25] . The 
union of all Farey classes up to a certain order J has the 
cardinality of the set of peaks for a spectrum of maxi- 
mum rank J. Now, the cardinality of Pj is know from 
number theory to be Euler's totient function, cf>(j) [2"o] . 
Therefore, the degeneracy D of the energy minima for a 
system with a maximum peak rank J is |26j 

J „ 

3 = 1 

which, from Eq. [2] scales as D ~ 2R/a. 

Finally, we recall j3[ [18] that the order j\, j% of the 
peaks bracketing a minimum relates to its structure in 
a straighforward way: the helix corresponds to a rhom- 
bic lattice where each particle has its nearest neighbors 
at axial displacements of iaji , ±aj2 and second nearest 
neighbors at ±a(ji or ±a(ji —32) [IB]. Also, ji and 
J2 give the number of crossing secondary spirals (paras- 
tichies) needed to cover the lattice by connecting nearest 
neighbors. 
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FIG. 3: Lattice energy V(il) versus divergence angle for suc- 
cessively halving values of a/R starting from 0.5 (using dipole 
dipole interaction, e = p 2 /a 2 , where p is the magnetic dipole). 
Notice the proliferation of peaks as a/R decreases. Repro- 
duced from |22j . 

For completeness, let us now follow Levitov [T8] . [23] . [24 ] . 
and consider the adiabatic evolution of our system as 
the linear density is increased. As new sets of maxima 
and minima emerge, the true minimum goes through 
a series of quasi-bifurcations, the consequence of an 
elusive symmetry whose explanation goes beyond our 
scope. Suffice it to say that the system evolves quasi- 
statically from one of these optimal fl to another as 
R/a increases, asymptoting to the golden angle fix = 
27r/ (r +1) [r = (l + \/5) /2] , ubiquitous in botany, as 
each minimum is bracketed by peaks whose ranks, be- 
cause of the Farey tree structure described above, are 
consecutive elements of the Fibonacci sequence. Occa- 
sional "wrong turns" at later stages, will not shift the 
convergence too far from the golden angle, yet the Fi- 
bonacci structure would be lost. However if one or two 
consecutive wrong turns happen at the second or second 
and third bifurcations the system will converge to the al- 
ternative angles of second or third phyllotaxis, given by 
Eq. 0. 

We have only surveyed so far spiraling lattices gener- 
ated by a single helix. A straightforward generalization 
gives multijugate phyllotaxis, when two or more elements 
grow at the same axial coordinate [IMS] ■ This case, which 
Levitov does not explore, can be easily mathematically 
reduced to monojugate case, by considering two or more 
replicas of the phyllotactic lattice as in Fig. [2] In our ex- 
perimental realization we restrict ourselves to the mono- 
jugate phyllotaxis, and explore multijugate only numeri- 
cally. 



III. A MAGNETIC CACTUS 

There is a long history of experimental reproductions 
of phyllotactic patterns. Recently, Doady et al. described 
phyllotaxis in terms of dynamical systems and then veri- 
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FIG. 4: Experimental apparatus. Top: Each unit of the mag- 
netic cactus consists of a magnet element and a unit ring 
secured to a central axis. The ring diameter d is 2.2 cm. 
Bottom: a schematic representation of the mounted magnetic 
cactus and surrounding measurement devices. The viewer's 
eye is restricted by the viewing slit and the reference wires. 
Measurements are taken directly from the dividing head. 



fied it experimentally by examining dynamical processes 
in droplets of ferro- fluid [57] . But even more than a cen- 
tury ago, Airy showed that phyllotaxis emerged in opti- 
mal packing of hard spheres connected by a rubber band, 
once the band was twisted to increase density [55] , 



Here we expand on what was announced in a recently 
published Letter |22) : we verify experimentally the as- 
sumptions of Levitov's energetic model, by studying the 
low energy configurations of interacting magnets stacked 
evenly-spaced and free to rotate around a common axis. 
We constructed a mechanical system that it is free to 
explore the three angles of botanical phyllotaxis (Eq. [IJ . 



A. Experimental Apparatus 

We built a magnetic cactus by mounting permanent 
magnets (spines) on stacked co-axial bearings (a stem) 
which are free to rotate about a central axis, as in Fig. [4] 
All the magnets point outward, to produce a repulsive 
interaction between all magnet pairs. To avoid effects 
of gravity, the apparatus rests in the vertical position, 
and is non-magnetic. We built two different versions, the 
second with magnets twice as long as the first, as to have 
a larger effective radius which gives three rather than two 
stable structures. 

We employed cylindrical permanent iron-neodymium 
magnets, 1.2 cm long and 0.6 cm in diameter. They are 
mounted on fifty aluminum rings of 2.2 cm outer diam- 
eter, each affixed to a non-magnetic radial ball bearing 
(acetal/silicon, Nordex) as in Fig.|4j These unit rings are 
evenly spaced on an aluminum rod in a stacked structure 
of 39.9 cm axial length. 

At static equilibrium we measure separation angle be- 
tween each magnet element, by rotating the cactus un- 
til a magnet element aligns with the reference wires. A 
telescope and a vertical viewing slit accompanied by two 
vertical reference wires assist in data acquisition. 



B. Annealing 

By measuring the dipole-dipole interaction between an 
individual magnet pair, we can reconstruct the curve of 
the lattice energy V(fi) as a function of the angular offset 
n between magnets. We find that the first arrangement, 
with short magnets, admits two minima, given by the 
angles of Eq. [T] for p = 1,2. The second arrangement, 
with long magnets, has three minima corresponding to 
the angles of Eq. [I]p = 1, 2, 3, one of which (p — 3) is a 
weak metastablc minimum. These divergence angles of 
stable helices are very close to those predicted by phyl- 
lotaxis, of Eq. [I] and are all accessible by experimental 
procedure described below. 

Before every data acquisition, the cactus is disordered 
and then athermally annealed into a low-energy state. 
The protocol involves repeatedly winding the bottom- 
most magnet to generate an ever-tightening spiral, un- 
til an explosive release of energy disorders the lattice. 
Next, an independent external magnet is oscillated in 
small circular motions near randomly chosen points while 
the cylinder as a whole is slowly rotated, to further ran- 
domize magnet orientations. After 10-30 second of me- 
chanical annealing through applied vibrations, the sys- 
tem consistently enters a robust ordered state which does 
not anneal further on experimental timescales. 



C. Results 

Figure [5] reports the experimental results for both ar- 
rangements by plotting the measured angle between con- 
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FIG. 5: Top: A 3-D rendering of the experimental data and 
the corresponding Bravais lattice of the magnetic cactus an- 
nealed in a spiral configuration of divergence angle Qi and 
parastichies (2,3) (blue and red dashed lines), Fibonacci num- 
bers. The bottom two panels show the experimentally mea- 
sured angular offsets Q between successive magnets for mag- 
netic cacti with short (middle) and long (bottom) magnets, 
plotted versus magnet index, which simply counts the num- 
ber of magnets along the axis. Flat regions are perfect spirals 
while steps are boundaries between different phyllotactic do- 
mains. The dotted lines give the phyllotactic angles fli, Q2, 
SI3 and 2n — SI2 defined in the text, whereas the dashed lines 
are minima of the magnetic lattice energy (insets) calculated 
by interpolating the measured pair-wise magnet-magnet in- 
teraction. Data reproduced from Ref. [22] , 
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FIG. 6: A kink between domains of first and second phyl- 
lotaxis. From top to bottom: a 3-D rendering, and its un- 
wound 2-D Bravais lattice from numerical simulations. Below, 
the same kink plotted as angle increments between successive 
magnets for the experimental system (crosses) and numerical 
simulation (circles). 



secutive magnets. The more narrow (short-magnet) cac- 
tus self-organizes into the spirals with divergence angles 
precisely reproducing those of first phyllotaxis, Q = fii, 
and second phyllotaxis, f2 = il 2 , as in Eq. 0. When the 
results are represented in a 2-D lattice, as in the top of 
Fig. [5] parastichies can be drawn. As parastichial num- 
bers for ft — fli we find the Fibonacci numbers (2,3), 
and for ft — fi 2 > the Lucas numbers (3,4), as seen also 
in botany. The larger-radius system also forms first and 
second Phyllotaxis helices, as well as limited domains of 
third phyllotaxis [with O = ^3 and Lucas numbers (1,4)], 
bracketed by domains of second phyllotaxis. The insets of 
Fig. [5] show the magnetic interaction energy V (f2) of the 
lattice obtained by interpolating measured values for the 
pair-wise magnet-magnet interaction, plotted as a func- 
tion of divergence angle ft. As we cans see, local minima 
correspond to phyllotactic angles. 

Figure [5] also shows that in many instances the system 
fragments into two or three distinct domains whose do- 
main walls always share a common parastichy, as seen in 
botany [UQ3], and as expected in physics for a quasi-one- 
dimensional degenerate system. We have computed nu- 
merically one such transition via dynamical simulations 
in a velocity- Verlet algorithm, in the following way: we 
start from a crude static step-like kink as an initial condi- 
tion, and allow it to radiate energy in the form of phonons 
waves until it stabilizes in a kink with superimposed vi- 
brations; we then average this configuration over time, 
to remove these residual oscillations. When the result is 
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FIG. 7: A numerically calculated kink in a system of dipoles 
of high degeneracy (seven minima, a/R — 0.15). The kink 
separates domains with parastichy numbers (4,5) and (5,6) 
and divergence angles of 1.38 and 1.13 radians. The top and 
middle panels give its three-dimensional rendering and angu- 
lar shift fi versus the axial magnet index. The two domains 
correspond to minima bracketed by peaks at Q/2tv = 1/4, 1/5 
and Q/2n = 1/5, 1/6 of the lattice energy, given in the bottom 
panel, where e = p 2 /a 2 , p being the magnetic dipole. 



used as new initial conditions, it proves to be a static 
kink. Fig. [6] reports our numerical results for a kink in a 
system whose size and interaction reproduces the physi- 
cal realization of the magnetic cactus, along with the ex- 
perimental data for such a kink. The match is essentially 
perfect, indicating that the dissipative (i.e. frictional) 
forces neglected in our model do not significantly affect 
the static configurations. We apply the same numerical 
procedure to calculate a kink in a system of larger degen- 
eracy, among domains which are absent in our physical 
realization. We use a smaller a/R ratio and a differ- 
ent interaction between particles (ideal dipole instead of 
physical dipole). The result shown in Fig. [7] reproduces 
the same qualitative shape of the previous, lower degen- 
eracy case. Similar kinks are present also in a fully un- 
constrained cactus, one in which the particles can move 
along the axis, and are found in early generations of our 
structural genetic algoritms (see below). Finally, these 
kinks can travel as novel topological solitons, with a rich 
phenomenology that is explained elsewhere [22J, [29] . 
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FIG. 8: Measured angular offsets between successive magnets 
for magnetic cacti with long magnets, plotted versus mag- 
net index, showing symmetric kink/anti-kink domain bound- 
aries. Dotted lines give the phyllotactic angles fii, Q.2, SI3 
and 2ir — SI2 defined in the text. Dashed lines are minima of 
the magnetic lattice energy (inset) calculated by interpolating 
the measured pair-wise magnet-magnet interaction. 



Finally, in the system with longer magnets, we occa- 
sionally found intriguing yet hard-to-interpret configura- 
tions that contain two nearly reflection symmetric do- 
main boundaries. Fig. [8] reports two such configurations, 
measured in independent experimental runs. Although 
we do not have a firm explanation for these structures, 
we speculate that they form as frozen-in soliton waves 
that initiated symmetrically at both ends of the struc- 
ture, upon release of the wound-up elastic energy during 
initial preparation. Indeed an analytical, continuum the- 
ory for phyllotactic solitons which we have developed re- 
cently, and which explains results of the dynamical symu- 
lations also supports the existence of similar frozen-in 
pulses [29] , 



IV. FULLY UNCONSTRAINED CACTUS: 
STRUCTURAL GENETIC ALGORITHM 

Our experimental apparatus is not fully unconstrained: 
the axial coordinates of the dipoles are fixed, and so 
only the azimuthal movement is allowed. While this is 
an huge improvement toward the original helical con- 
straint of Levitov, many (most) physical systems that 
could manifest phyllotactic patterns do not posses such 
a lesser azimuthal constraint. To corroborate and extend 
our experimental results to a completely unconstrained 
system, we seek the energy minimum in a set of repulsive 



7 



particles that can move axially as well as angularly on a 
cylindrical surface, via a non-local numerical optimiza- 
tion. To this purpose, we developed a structural genetic 
algorithm. 

A. Genetic Algorithm 

A genetic algorithm is a method of optimization that 
mimics evolution to find the absolute minimum in a func- 
tion which shows a large number of metastable minima. 
The coordinates of the energy functions are called genes, 
and a set of genes is a particular specification of value 
for those variables. The routine typically starts with a 
set of "parents" , or specific points in the domain of the 
energy function. At each stage of the routine, parents 
"mate" to produce children via exchange of genes: a sub- 
set of the coordinates of each of the two configurations 
are swapped, therefore generating new points in the en- 
ergy domain, called children. Each of those children is 
then locally relaxed to a minimum via a local search. 
The new population of parents and children undergoes 
genetic selection and only the fittest (the lowest energy 
ones) form a new population. 

There are many different implementations of this gen- 
eral idea: care is taken not to lose genetic diversity during 
selection, to avoid a population of almost identical repli- 
cas; that is usually achieved with a more or less skilled 
genetic selection, which might retain less geneticaly fit 
individuals, and often by introducing mutations in the 
form of random alteration of the gene sequence, which 
would opefully prevent the routine from gettting stuck 
around a metastable region. Choice of parameterization 
of the structures (genes) and mating (crossover) is crucial 
to the performance of the algorithm. 

About fifteen years ago, Deaven and Ho [30] introduced 
a so-called structural genetic algorithm, which proved 
particularly efficient in minimazing the energy of physical 
structures, as it allows for physical intuition in defining 
the genes and mating procedure. With it, they found the 
Ceo fullerene structure as a ground state of 60 carbon 
atoms interacting with suitable atomic potentials [30] 
and solved the celebrated Thomson problem of repulsive 
charges on a sphere [31], a task quite similar to ours. 
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FIG. 9: Relative energy of the fittest member of the popula- 
tion in every generation. Early generations return very high 
energy configurations corresponding to disordered metastable 
states. Interrmediate generations show populations of phyl- 
lotactic domains separated by kinks between. Finally, the 
algorithm converges to a single crystalline domain in the bulk 
(deformations at the boundaries simply accomodate the sys- 
tem to the confining potential). 



an external axial square-well of width L, which sets the 
length of the cylinder, and hence the density. The choice 
of the pairwise interaction is not fundamental, as long as 
it is long ranged, repulsive, and well behaved |22j : our 
particular choice simply speeds up the computation. 

We generate the first population randomly. At each 
step, we randomly couple mates, and excahge their genes 
employing the following mating procedure: we order the 
genes by increasing axial coordinates Z\ < z 2 < • • ■ < 
zioi and swap the first 1 < n < 101 genes, where n is 
a random number, between randomly selected parents. 
The children obtained in this way are then relaxed to a 
stable structure via a standard conjugate gradient algo- 
rithm. We then prepare the new generation by selecting 
the lowest energy individuals in the population of parents 
and children, yet making sure that the energy difference 
between members does not fall below a certain threshold, 
to preserve genetic diversity: when new children cannot 
produce a new population of 10 in accordance with the 
energy threshold, we introduce mutations by randomly 
altering a certain number of members. 



B. Our Algorithm 

In our implementation we use a population of 10 mem- 
bers. Each member reppresents a configuration of 101 
particles on the cylinder: more esplicitly, the genetic 
structure of each member P^, k = 1, . . . , 10, of the popu- 
lation is a set of variables, or = {6^, z.;}*£\ which spec- 
ifies, in cylindrical coordinates, the positions of the par- 
ticles composing its structure. The particles interact via 
a pair-wise inverse quadratic repulsion V = V (r /r) 2 , 
where r is the three dimensional distance between parti- 
cles; we introduce a confining potential in the form of 



C. Numerical Results 

During the structural evolution, the earliest popula- 
tions contains metastable disordered states. Members of 
intermediate populations show kinks between domains of 
different divergence angle, configurations which are also 
seen experimentally. After fifteen to twenty generations, 
the algorithm typically converges to a single crystalline 
domain. 

Figure [9] reports the energy of the fittest member of 
the population at each generation in a typical run, show- 
ing a punctuated-equilibrium evolution where the most- 
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FIG. 10: Numerical optimization via structural genetic al- 
gorithm for N = 101 repulsive particles [V = V (r /r) 2 ] 
constrained to a cylindrical surface of length L and radius 
R — 1.65L/7V. The resulting 2-D Bravais lattice has a nearly 
constant axial separation Az = Zi+i — Zi (top) and angular 
divergence fl between successive particles (bottom), neglect- 
ing fringe effects at the border of the potential well. In the 
bulk, particles self-organize on a single spiral of divergence 
= £7i . Oscillations at the boundaries are due to the effect 
of the confining potential. 



fit structure progressively decreases in energy in intermit- 
tent steps separated by plateaux. The final converged re- 
sults form well-defined two-dimensional cylindrical crys- 
tals away from boundaries. 



Figure 10 shows the crystalline structure to which the 
algorithm converges, for R = 1.Q5L/N: a single spi- 
ral with ft = Hi, as defined in Eqn. (IT]), correspond- 
ing to first phyllotaxis with parastichies (1,2). A plot of 
Az = Zi + \ — Zi returns the value L/N in the bulk, which 
implies a single generative spiral. This choice of RN/L 
corresponds to a density close to that of our experimental 
apparatus. 



V. MULTIJUGATE PHYLLOTAXIS 

For highly degenerate systems the genetic algorithm 
returns configurations with more than one generative 
spiral, corresponding to what in botany is called mul- 
tijugate phyllotaxis [2 . We have seen before that he- 
lices make cylindrically symmetric lattices. On the other 
hand, every cylindrically symmetric lattice can be de- 
composed into a suitable number of equispaced genera- 
tive spirals [2j [9] . That is accomplished by discretizing 
the cylinder along its axis into equally spaced rings and 
then assigning at each ring n sites, equally spaced and 
separated by a 2n/n angular shift. As before, each ring 
is shifted consecutively by a divergence angle il. The case 




FIG. 11: Top: 3-D schematics and 2-D lattice for a 2-jugate 
configuration. Bottom: Lattice energy versus angular offset 
for n-jugate configurations in a system of repulsive dipoles 
(e = p 2 /a 2 , where p is the magnetic dipole) for 1-jugate 
(black), 2-jugate (red), 3-jugate (green), and 4-jugate (blue). 



n = 2 is shown at the top of Fig. [TTJ The case n = 1 is 
shown in our experimental arrangement of Fig. [TJ 

By decomposing the n-jugate cylindrical lattice into 
n lateral replicas of single-spiral lattice, as in Fig. [2j the 
reader is easily convinced that multijugate phyllotaxis re- 
duces to the previously described monojugate case. All 
the considerations above apply, provided that one now 
takes the periodicity to be 2ir/n, and the distance be- 
tween rings to be na (with, as before, a = L/N). It 
follows that a n-jugate configuration will have local max- 
ima in 2m/ n, i = 1 . . . n and, following the discussion of 
section II, one finds that there will be other local maxima 
corresponding to angles 27r/n x i/j when j < [[J/n]], and 
J is the maximum rank given by Eq. [2j 

Note now that if two multijugate lattices of jugation 
ri, n' have a peak in the commensurate angle 2m/ j, then 
the energy of the peak is the same, as is shown in Fig. [TTJ 
bottom, which compares the plots of the energy of such 
an arrangement for different values of n. In fact both 
configurations correspond to particles facing each other 
after j/n and j/n' rings, and therefore at the same dis- 
tance na x j/n = n'a x j/n' = ja, independent of n oi n' . 
For small n the minima in the energy of n-jugate config- 
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urations essentially degenerate with the monojugate one 
previously explored. For large n they have higher en- 
ergy. If a/R is small enough, the threshold is n > J, as 
interaction between particles on the same ring become 
comparable to those facing in the minimal monojugate 
peak. 

VI. CONCLUSION 

We have studied the lowest energy configurations of 
repulsive particles on cylindrical surfaces, both experi- 
mentally and numerically. We have found that they cor- 
respond to the spiraling lattices seen in the phyllotaxis 
of living beings, both monojugate and multijugatc. By 



establishing experimentally and numerically that phyl- 
lotactic point lattices are ground states in the very gen- 
eral geometric scenario of unconstrained repulsive parti- 
cles on cylinders, we have opened the study of phyllotaxis 
to a much wider range of annealable physical systems 
where the particles could be electrons, adatoms, ions, 
dipolar molecules, nanoparticles, etc. constrained by ex- 
ternal potentials. 

Unlike plants, these multifarious, non-biological Phyl- 
lotactic systems could access various degrees of dynamics, 
providing new phenomenology well beyond that available 
to over-damped, adiabatic botany. We have reported 
elsewhere [22] on the dynamical richness of this physical 
phyllotaxis, including classical rotons and a large family 
of novel, inter-converting topological solitons. 
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